See https://psyarxiv.com/5ygtc/
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
df <- haven::read_sav("data/Data_Study 1.sav") |>
mutate_all(as.numeric) |>
filter(!is.na(exclude) & exclude == 0) |>
mutate(gen = as.character(ifelse(gen == 1, "Male", ifelse(gen == 2, "Female", "Other"))))
data <- df |>
select(matches("npi[[:digit:]]"),
matches("narq[[:digit:]]"),
matches("pni[[:digit:]]"),
matches("hn[[:digit:]]"),
matches("dt[[:digit:]]"),
matches("dsm[[:digit:]]")) |>
select(-ends_with("r")) |>
normalize(verbose=FALSE) |>
mutate(across(everything(), as.numeric))
# names(df)
# data
# summary(data)
paste0(
"Data from the [study 1](https://osf.io/gp6a4/) (Weidmann et al.), downloaded from OSF, included ",
report::report_participants(df, age = "age", gender = "gen", race = NA),
"."
)
[1] “Data from the study 1 (Weidmann et al.), downloaded from OSF, included 5736 participants (Mean age = 21.3, SD = 6.9, range: [18, 75], 2.6% missing; Gender: 64.1% women, 32.5% men, 0.40% non-binary, 3.07% missing).”
plot_hist <- function(data, x) {
data |>
select(starts_with(x)) |>
pivot_longer(everything()) |>
filter(!is.na(value)) |>
ggplot(aes(x = value)) +
geom_histogram(aes(fill=name), alpha = 0.3, position = "dodge") +
theme(legend.position = "none")
}
patchwork::wrap_plots(
plot_hist(data, "npi"),
plot_hist(data, "narq"),
plot_hist(data, "pni"),
plot_hist(data, "hn"),
plot_hist(data, "dt"),
plot_hist(data, "dsm")
)
plot_hist <- function(data, x) {
data |>
select(starts_with(x)) |>
pivot_longer(everything()) |>
filter(!is.na(value)) |>
ggplot(aes(x = value)) +
geom_histogram(aes(fill=name), alpha = 0.3, position = "dodge") +
theme(legend.position = "none")
}
patchwork::wrap_plots(
plot_hist(data, "npi"),
plot_hist(data, "narq"),
plot_hist(data, "pni"),
plot_hist(data, "hn"),
plot_hist(data, "dt"),
plot_hist(data, "dsm")
)
# r <- cor(data, use = "pairwise.complete.obs")
r <- correlation(data)
r |>
arrange(desc(abs(r))) |>
head()
## # Correlation Matrix (pearson-method)
##
## Parameter1 | Parameter2 | r | 95% CI | t | df | p
## ------------------------------------------------------------------------
## narq6 | narq9 | 0.70 | [0.69, 0.71] | 67.89 | 4803 | < .001***
## pni8 | pni16 | 0.70 | [0.68, 0.71] | 61.86 | 4021 | < .001***
## pni8 | pni40 | 0.65 | [0.63, 0.67] | 53.91 | 4027 | < .001***
## pni30 | pni36 | 0.65 | [0.63, 0.66] | 53.94 | 4032 | < .001***
## pni31 | pni45 | 0.65 | [0.63, 0.66] | 53.84 | 4022 | < .001***
## narq11 | pni29 | 0.64 | [0.62, 0.66] | 50.53 | 3625 | < .001***
##
## p-value adjustment method: Holm (1979)
## Observations: 3627-4805
n <- parameters::n_factors(data, n_max = 15)
n
## # Method Agreement Procedure:
##
## The choice of 1 dimensions is supported by 2 (28.57%) methods out of 7 (Acceleration factor, VSS complexity 1).
plot(n)
efa1 <- parameters::factor_analysis(data, n=1, sort = TRUE)
efa3 <- parameters::factor_analysis(data, n=3, rotation = "oblimin", sort = TRUE)
efa3_varimax <- parameters::factor_analysis(data, n=3, rotation = "varimax", sort = TRUE)
efa3_equamax <- parameters::factor_analysis(data, n=3, rotation = "equamax", sort = TRUE)
efa3_bentlerQ <- parameters::factor_analysis(data, n=3, rotation = "bentlerQ", sort = TRUE)
wrap_plots(plot(efa3), plot(efa3_varimax), plot(efa3_equamax), plot(efa3_bentlerQ))